The general harmonic-oscillator brackets: 
compact expression, symmetries, sums 
and Fortran code 



G.P. Kamuntavicius ^, R.K. Kalinauskas^, B.R. Barrett 
S. Mickevicius ^, and D. Germanas^ 

^ Vytautas Magnus University, Kaunas LT-3000 Lithuania 
^Institute of Physics, Vilnius LT-2600 Lithuania 
"^Department of Physics, University of Arizona, Tucson, Arizona 85721 



Abstract 

We present a very simple expression and a Fortran code for the fast and precise 
calculation of three-dimensional harmonic-oscillator transformation brackets. The 
complete system of symmetries for the brackets along with analytical expressions 
for sums, containing products of two and three brackets, is given. 
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1 Introduction 



A basis of liarmonic-oscillator functions has proven to be extremely useful and 
efficient in describing compact quantum systems, such as nucleons in atomic 
nuclei and quarks in hadrons. The traditional applications, such as the nu- 
clear shell-model, however, are based on a model Hamiltonian with individual 
one-particle variables. Hence, the corresponding model wave functions, which 
are dependent on one-particle coordinates, are not translationally invariant 
and cannot represent the wave function of such a system in proper way be- 
cause the center of mass (cm.) of a free nucleus (or a free hadron) must be 
described by an exponential function, corresponding to a freely-moving point 
mass. This need not to be a problem in the case when the expression for the 
realistic wave function in an harmonic-oscillator basis is known. In fact, hav- 
ing this expression, at least two possibilities exist to find a solution of the 
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problem of translational invariance of the wave function. The first approach is 
based on constructing the superpositions of shell-model states which possess 
a fixed state for the cm. of the system, see [1]. Utilizing a basis of this kind, 
we do not need to worry about the cm. state because all operators for ob- 
servables are translationally invariant so that the final result is independent 
of this state. The second possibility is based on a direct construction of the 
many-fermion wave function, which is independent of the cm. coordinate, see 
[2], [3], [4]. In this case, the harmonic-oscillator basis set, in terms of intrinsic 
(Jacobi) coordinates, is necessary. By a set of Jacobi coordinates for a system 
of N particles, we mean the N — 1 independent vectors that each represents 
the displacement of the cm. of two different subsystems. For N > 2, there 
exists more than one set of Jacobi coordinates that can be assigned to an 
A'"-particle system. The different possible sets of Jacobi coordinates can be 
identified by use of a Jacobi tree (see, for example, [5]) and are related to 
each other by orthogonal transformations. In general, when the transforma- 
tion from one set of Jacobi coordinates to another is performed, one obtains an 
expansion for the wave function containing an infinite number of terms. Only 
the set of harmonic-oscillator functions can be chosen in such a manner that 
the transformation from one set of Jacobi coordinates to another results in a 
corresponding expansion with a finite number of terms. In any of the above 
mentioned approaches, the essential feature is the Talmi-Moshinsky transfor- 
mation [6], [7], [8] and corresponding harmonic-oscillator brackets (HOBs). 
Historically, the first application of this procedure was transformation of the 
product of two harmonic-oscillator functions from single-particle coordinates 
to relative motion and cm. coordinates, a reduction that has proven to be 
very useful for the evaluation of two-body matrix elements. Because the HOBs 
are constantly employed in various model calculations of nuclear and hadrons 
structure, it is desirable to have a simple and efficient method for calculating 
them. Many papers have been devoted to the study of the Talmi-Moshinsky 
transformation and brackets, and various methods for the calculation of these 
brackets and several explicit expressions for them are described in the litera- 
ture (see [9] and references therein, [10] and references therein). In this paper 
we present a complete system of symmetries for these brackets, a very simple 
expression for the HOBs based on the result of [10], new expressions for the 
sums of products of HOBs and a computer code, written in Fortran, which 
calculates the HOBs quickly and precisely. 
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2 The definition 



First, let us consider the HOBs, defined in the following way: 

\eih (ri) , 62/2 (r2) : AA) = ^ {EL, el : A \eih, 62/2 : A)^ \EL (R) , el (r) : AA) , 

EL,el 

(1) 

or, in other words: 



{EL, el : A |ei/i,e2/2 : A)^ 

= ^aVi ^ / / ^^"^"^ ' ^""^ ■ "^^1 ""'^^ ^""^^ ' "^2^2 (""^^ ■ ■ 

A 

Here c? is a nonnegative real number (sec Eq. (7)). The two-particle wave 
functions with bound momenta are defined as: 



\eili (ri) , 62/2 (r2) : AA) = {(p^.h (ri) ® 'Pe2i2 i'^2)}Ax 

= {hmi,l2m2 |AA) 0eiiimi (ri) 0e2i2m2 (r2)(3) 



mim2 



and the properly orthonormalized harmonic-oscillator function is given by 

exp (-rV2) r^4'+V2) (^2j ^^/^^ ^ 



^i'dm (r) = (-1)' 



2 (n!) 
r(n + Z + 3/2) 



(4) 

where the corresponding dimensionless eigenvalue equals (e -1-3/2), n — 
(e — 0/2 = 0, 1,2, .... We prefer the quantum number e = 2n -\- I, rather 
than n, due to conservation of the total oscillator energy on both sides of the 
bracket: 

ei + e2 = £; + e. (5) 

Obviously, this relation gives the requirement for the angular momenta: 

(-l)'i+'2 = (-1)^+' . (6) 

Let us next define the transformation of variables present in the expression 
for the brackets in the following way: 



/ 
V 



R 

r 



n~ \ ( V \ 

V r2 J 



1 

l+d 



l+d / 



\J \+d \l l+d 



V Vl+<i J V 

(7) 



R 

r 
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Wc prefer the above order for the variables and the couphng of the angular mo- 
menta because this produces an orthogonal and, at the same time, symmetric 
transformation matrix, hence, giving higher symmetry to the brackets: 



{EL, el : A |eiZi, 62^2 : A)^ = (ei/i, 62/2 : A \EL, el : A) 



d ■ 



(8) 



Any orthogonal matrix of second order can be represented in form (7). If 
problems arise, the matrix can be easily rewritten in this form. The complete 
solution of this problem is as follows: In general, the matrix in (7) has the 



form 



^ sin 6* cos 6' 



cos 9 — sin ^ 



and possesses the required distribution of signs, i.e.. 



+ + 



, only in the case < 9 < 7i/2. When 7t/2 < 9 < n, the signs change 



so that the matrix looks like 



, but by an elementary transformation 



V 



of variables, namely R — > — r and r — > R, we can transform it to the previous 
form, i.e., only one minus sign in the lower right corner. If tt < 9 < 37r/2, 

/_ _\ 

distribution of signs is and the necessary transformation of variables 

\-+J 

is R — > — R and r — > — r. The last case, when 37r/2 < 9 < 27t, gives the matrix 

l-A 

; the correcting transformation is now R — > r and r — > — R. When 

V++; 

the transformation is given by a nonsymmetric orthogonal matrix, the HOBs 
must, in general, be written in the form 



EL el 
R r 



A 



eih 62/2 



ri r2 



A 



(9) 



which clearly demonstrates the correspondence between quantum numbers 
and variables. Some authors define these brackets without referencing this 
correspondence; in such a case, one must be careful, when employing some of 
their symmetries. 



The HOBs, defined as above, are the real entries of an orthogonal matrix. 
Consequently, they obey usual orthogonality conditions: 



E 

EL,el 



{eih, 62/2 : A \EL, el : A)^ {EL, el : A \e^l^, e^l^ : A)^ = ^eiZi.eiii^eaia.e^i^- 



(10) 
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^ {EL, el : A |eiZi, 62^2 : A)<i (ei^i, 62^2 : A \E' L\ e'l' : A)^ = SEL,E'L'Sei,e'i'- 
61(1,62^2 

(11) 

Here and below we apply the compact expressions for the Kroneckcr deltas, 
e.g., Sel,e'L' = 5e,e'Sl,l'- The symmetry properties of the coefficients are: 



(61/1,62/2 : A \EL,el : A)^ 

= {EL, el : A |eiZi,e2/2 : A)^ 

= (-1)^+'^ (e2/2,6iZi:A|eZ,£;L:A), 

= (-1)^-^ (62^2, 61Z1 : A \EL, el : A)^^ 

= (-1)^-'^ (61/1,62/2: A |eZ,£;L:A),/^. 



(12) 
(13) 
(14) 
(15) 



The symmetry defined in (12) follows from the definition of the coefficients 
and the requirement for the angular momenta (6); it is already given in (8). 
The next symmetry, (13), is not so trivial. It can be derived using the fact 
that the same transformation matrix connects not only the original, but also 
the transformed variables: 




V 



d 
1+d 

1~ 
1+d 



1 

1+d 



d 
1+d 



I 



\ 



r2 



(16) 



-ri 



Using this condition for corresponding wave functions, one immediately ob- 
tains the symmetry relation (13). The symmetries (14) and (15) are even more 
complicated. The first one, (14), is based on observation that 



V 



R 



/ 



V 



1 

1+d 



d 
1+d 



1 

1+d 



i V 
\ ^1 



(17) 



where the transformation matrix corresponds to the parameter value 1/d in- 
stead of the original value d. Let us ilhistratc the derivation of the symmetry 
in this case. The expression for the brackets can be represented in the form: 

(61/1,62/2 : A \EL,el : A)^ 

= ^aVi ^ ^^^^^ (""^^ ' "^2^2 ^""^^ • ' ^""^ '■ 



2A + 



J ^"^l^^^ {0ei/i (ri) 0e2/2 {^2)}tx {(pEL (R) (pel (r)} 



AA ■ 

1^ 



Now it is enough to reorder the variables according to both sides of (17). To 
do this, the following simple, well-known expressions are necessary: 



0eim (-r) = (-1)' 0eim (r) 



(19) 
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and 



{0eiii (ri) 062^2 (1^2) } 



AA 



[-1 



{062/2 (r2) ® (ri)} 



AA 



(20) 



Applying these expressions, the right-hand side of (18) can be rewritten as 
(61/1,62^2 : A \EL,el : A)^ 

= (-1)'^+'^-^+^ E / / dvidr, {cj>,,i, (r2) ® Ki, {v,)}l^ {cj>EL (R) ® <Pei (- 



= (-1)''-'' (62/2, eiZi : A \EL, el : A)y^ . 
The last symmetry, (15), follows from the transformation 



>2^2,ei/i : A \EL,el : A)^^^ 



(21) 




ri 



Vt5 ~VtE / V 



(22) 



The self-consistence of the transformations and, hence, the symmetries can 
be checked by applying two transformations, say (14) and (15), for the same 
bracket. The result yields an identity, as expected: 

(ei/i, 62/2 : A \EL, el : A)^ = (-1)'^""^ {e^k, e^h ■ A \EL, el : A),/^ 

= (-1)^-^ (-1)^-'^ (62/2, ei/i : A \el, EL : A), 
= (61/1,62/2 : A |£;L,e/ : A)^. (23) 

Hence, our definition of brackets is based on an assumption that variables in 
brackets are arranged in a fixed way, and the corresponding transformation 
matrix is given as in (7): 



(R,r|ri,r2),- \ _ nj: 

V 1+d V 



(24) 



3 Compact expression 



The HOBs have been considered previously by a number of authors; however, 
these results were quite complicated and led to expressions, whose structures 
are not transparent. In our opinion, the simplest known expression for the 
general oscillator bracket is the one derived by B.Buck and A. C. Merchant in 
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Ref.[10]: 



{EL, el : A |ei/i,e2/2 : A), 



X ^(m)! (ris)! {N)\ {n)\ [2 (m + Zi) + 1]!! [2 (712 + h) + 1]!! [2 (TV + L) + 1]!! [2 (n + + 1]! 



abcdlalb^c^d 



X 



[(2/, + 1) (2/fe + 1) (2/e + 1) (2/rf + 1: 



a\b\dd\ [2 (a + /„) + 1]!! [2 (6 + k) + 1]!! [2 (c + Q + 1]!! [2 {d + Q + 1]! 

la lb h 

x< Icldk) {Ia0lc0\ LO) {hOldOl 10) {Ia0h0\ hO) {Ic0ld0\ I2O) , (25) 
L I A 



where N ^ {E - L) /2, n = {e - I) /2, rii = (ci - h) /2, and 712 = 
(62 — ^2) /2. This expression for the HOBs is derived using harmonic-oscillator 
wave functions without the phase multiplier (—1)" present in our Eq. (4). The 
introduction of this modification results in a slightly different phase in the 
expression for HOBs. This phase equals (^—xj^+^+^i^^^ ^ ^_-^-^(L+i-h-i2)/2 ^ 

Although the above summation runs over eight indices, the real summation is 
only over five of them due to three independent constraints. This constrained 

summation is similar to other, known, modern expressions for the HOBs. For- 
mula (25) is very symmetric, and due to this symmetry, can be rewritten in 
the following way: 



{EL, el : A \eih,e2l2 ■ A)^ 

la h h 
Ic Id I2 

L I A 

X G {eih; eJa, ebk) G (62/2; eJc, eJd) G {EL; eJa, eJc) G {el; e^k, eJd) , (26) 



where 



G {eih; eJa, ebk) 



X 



(2/, + l)(24 + l)(/,0/50| /lO) 

ei-h \ I ei + h + 1 
ea - la] eb-lb \ea + la + i;eb + lb + i 



1/2 



(27) 
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and the coefficients 



/ 




(28) 



\na;nb 



are trinomials defined for tfie parameter values rii > na,nii . Moreover, the 
three parameters of these coefficients are all even or all odd at the same time. 
It should be noted that the well-known binomial coefficients can be expanded 
in trinomials, e.g.. 



In expression (26) we again have eight summation indices connected in four 
pairs, Cala (ct — ct; b, c, d), with the well-known relation between the oscillator 
energy and the angular momentum: is a nonnegative integer and 1^ — 
Cq, — 2, 1 or 0. The constraints amongst the energies are: 

Ca + Cc = E, eb + Cd^ e, 

Ca + 66 = 61, Be + 6^ = 62. (30) 

Due to relation (5), (£"-1-6 = 61 + 62), only three of them are independent. 
The best choice for the independent summation indices is 6^ and the four 
angular momenta la (a = a, b, c, d) . 



4 Sums of products 

The sums of HOBs occur when antisymmetrizing the translationally invariant 

wave function, i.e., the function of Jacobi and intrinsic (spin, isospin, etc.) 
variables, because any permutation of one-particle coordinates results in or- 
thogonal transformations of the Jacobian coordinates. These permutations 
result in different orthogonal transformations; therefore, according to (7), in 
transformations with different values of the parameters d. The keys to the 
derivation of the expressions for the sums of products are Eq. (18) and the 
possibility to represent any orthogonal transformation in the form (24). To 
illustrate this procedure, let us take two expressions for the HOBs, as given 



by Eq. (1), 

{(f>EL (R) «) (f>el (r)}^;, = J2 (^^' ■ ^ l^l^l' ^2^2 : A)d {0ei/i (ri) (8) 0e2/2 (r2)}AA ' 




(29) 



eih,e2h 



(31) 
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with different variables and different ortliogonal transformations, and multiply 
them, the left side to the left side and the right side to the right side. Next let 
us sum over the projection of the angular momentum and integrate over both 
coordinates to obtain: 

A 

J2 {EL, el : A \e,h, 62/2 : A)^ {E'L', e'l' : A \e[l[, e',l', : A),, 

eilie2l2e'^l'j^e'2l'2 

(32) 

Obviously, the variables R, r and ri,r2 are connected by the orthogonal 
transformation (R, r |ri, r2)^, Eq. (24), and variables R',r' and r^^r^ are 
connected by another orthogonal transformation (R',r' jr'j^jr'g)^, of the same 
structure. In the case when ri,r2 and r[,r2 are also connected by some or- 
thogonal transformation, say (ri,r2 \r[,r'2)^^^, both integrals on the left and 
the right sides of Eq. (32) are expressible as HOBs. Thus, Eq. (32) can be 
rewritten as 



{EL, el : A \E'L',e'l' : A) 



D 



J2 {EL, el : A \eili, 63/2 : A)^ {eih, 62/2 : A \e[l[, 62^2 : A) 



do 



X {e[l[,e'2l2 : A \E'L',e'l' : A)^, 



(33) 



The only remaining problem is the definition of the parameter D. Having in 
mind that all transformations of variables are orthogonal and well-defined, we 
can expand 



d 

l+d 
1 

l+d 




l+do 



dp 
l+do 




l+d' 



l+d' 

(34) 



Multiplying the above matrices, one obtains the form: 




D 
1+D 




(35) 



where 



'dd' 



1 -Vdd/ + Vd^(yd + y/d/)' 



(36) 
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and restrictions for values of the parameter do are as follows: If dd' < 1, then 
^/d^< (^Vd+Vd') / (l -Vdd'y,iidd' > l,thenv/4> (v^-l) / [Vd + Vd/). 
The sum in Eq. (34) of the product of three HOBs with different values of 
the parameter d can be simplified to sums of products of only two HOBs. 
To do this, the values of the brackets corresponding to the extreme values of 
the parameter d are necessary, i.e., = or — > oo. The transformations of 
variables in these cases are trivial and HOBs can be calculated directly from 
definition (18), yielding: 

{EL, el : A \eili,e2l2 ■ ^)d=0 = {-'i-)^'^^"^ ^eih,el3e2l2,EL, 

{EL, el : A \eih, 63^2 : A)d^oo = ^eih,ELSe2i2,ei- (37) 

Using the HOB with the parameter do ^ 00, one obtains the following ex- 
pressions (the value do ^ 00 is consistent only with dd' > 1) : 

{EL,el:A\E'L',e'l':A)^^^_^^^^^^^^^y 

= ^ (-1)'^ {EL, el : A \e,h, e^k : A), {e,h, e^k : A \E'L', e'l' : A),, . 

eili,e2l2 

(38) 

Using the HOB with the parameter do — (consistent only with dd' < 1), one 
obtains the result: 

{EL,el:K\E'L',e'l' -.A),, rjr\,(. rnrw^ 

= (-1)^'+^ (-1)'' (^^' : A |ei/i, 62/2 : A), {e,h, e^k : A \e'l', E'L' : A),, . 

(39) 

These sums appear, when calculating the matrix element of the permutation 
operator of the particle coordinates in the translationally invariant basis of 
oscillator functions, see [4]. They are consistent: at d' = 1/d in both cases one 
can easily obtain the normalization condition for HOBs. 



5 Fortran code 

Using the above results, we have also developed a Fortran code for the fast 
and precise calculation of the HOBs in large quantities. Such a computer code 
is needed because any nuclear calculation, including mentioned above anti- 
symmetrization of a translationally invariant basis, requires large matrices. 
As one can see from the expression for HOBs, Eqs. (26) - (28), the main 
elements of these expressions are the Clebsch-Gordan coefficients with zero 
angular momenta projections, the 9 — j symbols and trinomial coefficients. 
Our code is based on the observation that all group-theoretical expressions 
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can be represented as products or sums of products of binomial coefficients. 
Binomial coefficients arc far more acceptable for large calculations than are 
factorials. For example, to represent 50! exactly requires a mantissa with 53 
significant numbers. At the same time, the mantissa of the binomial re- 
quires only 15 significant numbers, hence it can be stored using real numbers 
of double precision. Having this in mind, for Clebsch-Gordan coefficients we 
use the expression: 



' 21 \ h + h + l + V 
M-k + l [ 21 + 1 , 



.-1/2 



I \ ({h + I2 + I) /2' 

Xh-h + l)/2)\ I , 

(40) 



For the 9 — j symbols we first express them in terms of 6 — j symbols, using 
the well-known formula: 



3l 32 J3 

h h h 
ki k2 k^ 



E(-ir(2a;+l) 



3i 32 is 

< 


> < 




h 


> < 


ki k2 


ks^ 

> 


^3 kj, X 

V / 




32 X 


k2 




_ X jl 


h 



(41) 



and then utilize the Q — j expression in terms of binomials: 



a b e 
dcf 



■.a+d+e+f 



a+c- f+1 



(a + b + e + l)(c + d + e + l] 



(a + c + f + l) (b + d + f + l) 



a+c+f 
2a 



X 



xE(-i) 

2 



) {a+c-f) ( 



c+d+eX / 2d \/a+e+b 
2d )[c+d~e)[ 2e 



( V 

\a+e-bj 



1/2 



1/2 



fb+f+d\f 2/ 
V 2/ J\b+f-d, 
fa+f-c\ ( f+c~a 

\ z J U+f-b-z 



/ a+c+f \ ^b+c+e~f+l+z\ • 
\a+d+f-e-zj V C+0-/+I 



(42) 



This formula for the Q—j symbol is well-balanced because every unit summed 
over has an equal number of binomials in the numerator and in the denomi- 
nator. This results in the highest precision for the sum. Our code starts the 
calculations of the HOBs by filling the arrays of the binomial and trinomial 
coefficients. For the binomials we use recurrence formulas with corrections tak- 
ing into account the fact that the binomial coefficients are integer numbers. 
For the three-dimensional array of trinomials we apply a completely analo- 
gous method. When both arrays are filled, the calculation of the HOBs is an 
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extremely fast and precise operation. In Table 1 we give values of 



Eo 




J2 {EL, el : A \eih, 62^2 : A)^ {eih, 62^2 : A \E'L', e'l' : A)^ 

61^1,62^2 

(43) 



which characterize the precision of our calculations. Equation (43) represents 
the sum of the absolute values of the deviations of the calculated normaliza- 
tion conditions for the HOBs from the exact values given by Eq. (11). The 
first sum runs over all free parameters of the normalization condition, i.e., 
E, L, e, /, E', L', e', I' and A, taking all allowed values between and £"0. The 
total number of HOBs calculated for a given value of Eq equals (Eq) . The 
processor time (in seconds), necessary to perform these calculations on a per- 
sonal computer, is listed under T {Eq) . 

Table 1 



Eq 


N{Eq) 


5(£o) 


T{Eq) 


3 


4,672 


0.2956* IQ-^* 


0.28 


4 


32, 054 


0.6832 * 10-13 


1.65 


5 


157.618 


0.5028 * IQ--'-- 


11.9 


6 


658, 000 


0.2394* 10-11 


74.7 


7 


2,298,144 


0.7336* 10-1° 


385 


8 


7, 165, 706 


0.7032 * 10-9 


1,760 


9 


19, 973, 174 


0.2057* 10-^ 


6,830 


10 


51,349,192 


0.1306* 10-^ 


25, 100 


11 


122,193,968 


0.2131* 10-5 


79, 900 


12 


273, 872, 766 


0.3481 * 10-^ 


248, 000 



As Table 1 clearly demonstrates, our method for calculating the HOBs is 
fast and universal and, hence, applicable to any calculations using a basis of 
many-particle harmonic-oscillator functions. Our Fortran code is available for 
general use and distribution, as described in the Conclusions. 



6 Conclusions 



Starting with the result of Ref.[10] for the General Harmonic-Oscillator Brack- 
ets (HOBs), we have simplified this earlier expression into a more compact and 
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highly symmetrical relationship, given in Eqs. (26)-(27). Our new result for the 
HOBs is well-suited for fast and precise numerical calculations of large quanti- 
ties of HOBs. The numerical procedure is based upon independent calculations 
of the Clebsch-Gordan coefficients of zero angular momentum projections, the 
6 — j symbols and the 9 — j symbols using binomial and "trinomial" coeffi- 
cients, which are precalculated and stored in appropriate arrays. The Fortran 
code, which we have written for computing the HOBs, is available for gen- 
eral use and can be obtained from the Vytautas Magnus University, Kaunas 
Lithuania, web server at |http: / / www.nuclear.physics.vdu.lt7. 



References 



[1] Neudachin V.G., and Smirnov Ju.F., Nucleonic associations in light atomic 
nuclei, Nauka, Moscow 1969, 414 p. (in Russian). 

[2] Vanagas V., Algebraic methods in nuclear theory, Mintis, Vilnius, 1971, 378 p. 
(in Russian). 

[3] Kamuntavicius G.P., Sov. Journ. Part, and Nuclei, 20, N2, 261 (1989). 

[4] Navratil R, Kamuntavicius G.P., Barrett B.R., Rhys. Rev. C, 61, 044001 (2000). 

[5] Smirnov Yu.F., Shitikova K.V., Sov. Journ. Part, and Nuclei, 8, 847 (1977). 

[6] Talmi I., Helv. Rhys. Acta, 25, 185 (1952). 

[7] Moshinsky M., Nucl. Rhys., 13, 104 (1959). 

[8] Smirnov Yu. F., Nucl. Phys., 27, 177 (1961); 39, 346 (1962). 

[9] Trlifai L., Rhys. Rev. C, 5, 1534 (1972). 

[10] Buck B., Merchant A.C., Nucl. Rhys. A600, 387 (1996). 



13 



